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Abstract 

A mathematical model for the quantitative analysis of cancer immune interaction, 
considering the role of antibodies has been proposed in this paper. The model is 
based on the clinical evidence, which states that antibodies can directly kill can- 
cerous cells [1]. The existence of transcritical and saddle-node bifurcation, which 
has been proved using Sotomayor theorem, provides strong biological implications. 
Through numerical simulations, it has been illustrated that under certain therapy 
(like monoclonal antibody therapy), which is capable of altering the parameters of 
the system, cancer-free state can be obtained. 

Key words: Cancerous cells, B-cells, Plasma cells, Antibodies, Transcritical 
bifurcation, Saddle-node bifurcation 



* Corresponding Author. 

Email addresses: yoonanet@gmail.com (Shiferaw Feyissa), 
sandofmaOiitr . ernet . in (Sandip Banerjee). 



Preprint submitted to Mathematical Biosciences 



1 Introduction 



Cancer is a leading cause of death worldwide. It is predicted to remain so for 
years to come unfortunately. It is a disease which is caused by the abnormal 
function of our own cells. Cancer is not a single disease; it comprises more 
than 200 diseases [2,3] which share common characteristics, that is, they are 
abnormal cells where the normal processes which regulate normal cell prolif- 
eration, differentiation and death (cell apoptosis) are interrupted. The causes 
of cancer are changes that cause normal cells to acquire abnormal functions. 
These causes may be the result of inherited mutation or environmental factors 
such as tobacco products, ultraviolet radiation, X-rays, chemicals, etc [3]. A 
normal cell can be transformed into a cancerous cell when certain genes are ac- 
tivated or inactivated because of these mutation and environmental factors. In 
cancerous cells the normal control system that prevent abnormal cell growth 
and differentiation and the invasion of other tissues and organs are disabled. 

The immune system, which is a complex network of cells, cytokines, lym- 
phoid tissues and organs that work together, helps the host in fighting against 
pathogenic micro organisms and cancerous cells. When cancerous cells prolif- 
erate to a detectable threshold number in a given physiological space of the 
human anatomy, the body's own immune system is triggered into a search 
and destroy mode [4]. The various components of the immune system interact 
with each other and the cancerous cells to prevent our body from cancer and 
they also help in preventing the occurrence, development and recurrence of 
cancer. But sometimes cancerous cells are able to overcome the limitations 
imposed by the immune system and are able to successfully proliferate and 
propagate in the human body. The defense mechanisms of the human body 
against external invaders and other pathogens has two major components; 
innate (or natural) and adaptive (or acquired) immune systems. The innate 
immune system comprises cells and mechanisms that defend the human body 
from infection by other organisms and pathogens in a non specific manner. 
Innate immune responses are stimulated by structures that are common to 
groups of related pathogenic agents and may not distinguish fine differences 
among foreign substances and it does not confer long lasting or protective 
immunity in our body. 

The major components of innate immune system are natural killer (NK) cells, 
dendritic and mast cells, macrophages and natural antibody producing cells. 
This immune system which is naturally present in our body even in the absence 
of external invaders or cancerous cells, provide immediate defense against in- 
fection. The biological mechanisms of innate immune system are integrated 
with that of adaptive immune system by stimulating and influencing the na- 
ture of the adaptive immune responses. The acquired immune responses are 
highly specific for a particular pathogen. These immune responses are charac- 
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terized by an exquisite specificity for distinct external invaders and pathogens 
and the ability to remember and respond more vigorously to repeated ex- 
posures to the same pathogens. Adaptive immune system responses against 
external invaders and pathogens are mediated by a special group of immune 
cells called lymphocytes through either the anti-body (humoral) mediated or 
cell-mediated immune responses [2,4]. The most important cells in acquired 
immune system are lymphocytes, which are small round cells that are found 
in blood, lymph and connective tissue. The two major components of lympho- 
cytes are the B-cells and the T-cells [5]. Both B-cells and T-cells are derived 
from the stem cells residing in the bone marrow. Cells which are destined 
to become T-lymphocyte undergo further differentiation in the thymus while 
the precursors of B-lymphocyte differentiate in other lymphoid organs. The 
B-lymphocytes are responsible for antibody mediated immune response and 
the T-lymphocytes are responsible for cell-mediated immune responses. 

The theoretical study of cancer-immune dynamics, which has a long history, 
has been done by many authors [6,7,8,9,10,11,12,13,14,15,16,17,18] . We first 
look through some of the existing mathematical models of the cancer-immune 
system interactions. Kuznetsov and Taylor [6] presented a mathematical model 
of the cytotoxic T-lymphocyte and responses to the growth of an immunogenic 
tumor. They studied the immune-stimulation of tumor growth, "sneaking 
through" of the tumor and the formation of a tumor dormant state. Through 
mathematical modeling Kirschner and Panneta [9] have illustrated the dy- 
namics between tumor cells, immune effecter cells and interleukin-2 (IL-2). 
Their efforts explain both short -term tumor oscillations in tumor size as well 
as long-term tumor relapses. Bodnar and U.Fory's [10] studied the periodic 
dynamics in the mathematical model of the immune system. In Murchuk's 
model of immune system dynamics, U.Fory's [11] presented the model of a 
general immune system reaction. The qualitative behavior of the solution to 
the model (and its application), along with many illustration of the recov- 
ery process, oscillations or lethal outcomes of the disease has been discussed. 
In [14], the authors expressed the spontaneous regression and progression of 
a malignant tumor system as a prey-predator like system, where the tumor 
(cancerous cells) is treated as prey and the cytotoxic T-lymphocytes as preda- 
tor. The deterministic model is extended to a stochastic one, allowing random 
fluctuations around the positive interior equilibrium. The stochastic stability 
properties of the model are investigated both analytically and numerically. 
Mallet and De Pillis [15] have presented a hybrid cellular automata partial 
differential equation model of moderate complexity to describe the interaction 
between tumor and the immune systems of the host. Chaplain et. al. [16] have 
explained the effect of time and space in tumor immunology using mathemat- 
ical model, that is, the spatio-temporal phenomena. The role of interleukin-2 
(IL-2) in tumor dynamics is illustrated through mathematical modeling (a 
modified version of the Kirshner-Panetta model) in [17], where the author has 
shown that interleukin-2 alone can cause the tumor cell population to regress. 
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In [18], the global dynamics of Kirshner-Panetta model is explored and under 
what conditions tumor clearance can be achieved, is obtained. 

In most of these tumor-immune system interaction models, the immune sys- 
tem comprises both antibody mediated and cell mediated immune responses. 
In such situations, it is difficult to identify which immune system plays greater 
role in the elimination or control of the cancerous cells. In this paper we have 
investigated the role of antibodies in the eradication of the cancerous cells 
using system of non linear differential equations. The motivation came from 
the fact that researchers at the University of Manchester along with their 
collaborators at the University of Southampton investigated how antibody 
treatments make cancerous cells kill themselves and found a previously undis- 
covered mechanism that could, in future, be even more effective in causing 
their death. It is known that when antibodies bind to cells, including cancer- 
ous cells, they can mark those targets for destruction by the body's immune 
system but Tim Illidge et. al. [1] have shown in their latest study that antibod- 
ies can kill cancerous cells directly. When the antibodies binds with cancerous 
cells, it causes lysomes (small acid containing sacs) inside the cell to swell and 
burst rapidly releasing their toxic contents with fatal results for cancerous 
cells, which is non-apoptotic in nature. 

Not much work has been done on the role of antibodies to eradicate cancerous 
cells through mathematical modeling. This may be due to the fact that no 
clinical evidence was available to support the fact that antibodies are actu- 
ally capable of killing cancerous cells directly. To the best of our knowledge, 
the role of antibodies to eradicate cancer by mathematically modeling the 
scenario were done by Dillman and Koziol [19], Kolev [20] and Dubey et. al. 
[21]. Dillman and Kziol proposed and developed a pharmacokinetic model for 
the quantitative analysis of dose-timecell survival curves devolving from in- 
fusions of the murine monoclonal antibody T101 into patients with chronic 
lymphocytic leukemia (CLL) and cutaneous T cell lymphoma (CTCL). Kolev 
has described a model of cellular tumor dynamics in competition with the 
immune system with the help of integrodifferential equations, where the role 
of antibodies has been taken into account. Dubey et. al. have considered many 
components of acquired immune response, namely, T Helper cells, Cytotoxic 
T-Cells, B-cells which secretes antibodies and their interaction with avascular 
cancer cells. They observed that under appropriate conditions this interaction 
is capable of controlling the growth of cancerous cells. 

In section 2, the whole biological scenario supported by schematic diagram 
is explained, followed by mathematical formulation of the model. Estimation 
of the system parameters are discussed in section 3 . Linear stability analy- 
sis, global stability analysis as well as bifurcation analysis of the system are 
explored in section 4. Section 5 deals with the numerical results and their 
biological implications. The paper ends with a conclusion. 
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2 Model Formulation 



The host immune system has the ability to produce some significant anti 
cancerous immune responses, one of which is the B-cells. The human body 
produces millions of different B-cells each day that circulates in the blood and 
lymphatic system performing the role of immune surveillance. The principal 
function of B-cells is to secrete antibodies against antigens. Prior to stimu- 
lation by either antigen or mitogen, B-cells are morphologically small cells. 
The binding of antigen to antigen receptors (i.e. antibodies) on B-cells can 
result in the activation and differentiation of small B-cells into large B-cells 
which secrete antibodies at a lower rate. A set of immunoglobulin molecules 
is present on the surface of unstimulated B-cells [22] . By binding to these im- 
munoglobulin receptors and with a possible second signal from an accessory 
cell such as the T-cell, the antigen stimulates the B-cell to divide and mature 
into terminal (non-dividing) antibody secreting cells called plasma cells. The 
plasma cells are most active in secreting antibodies at a much faster rate but 
large B-lymphocytes, which proliferate rapidly, also secrete antibody, albeit at 
a lower rate. Some of the large B-cells eventually revert back to small B-cells 
where they probably function as memory cells, which can respond vigorously 
to subsequent antigenic challenges [5]. The secreted antibodies then circulate 
in the blood and lymphatic system, and bind to the original antigen, marking 
them for elimination by several mechanisms, including activation of the com- 
plement system, promotion of phagocytosis via opsonization and mediation 
of antibody dependent cell-mediated cytotoxicity (ADDC) with effecter cells 
such as macrophage, NK cells and neutrophils. 

Based on the above biological scenario and from the clinical evidence that 
antibodies can directly kill cancerous cells [1], we present a schematic diagram 
(fig.l) to describe the interaction between the cancerous cells, large B-cells, 
plasma cells and antibodies. We now propose a mathematical model to de- 
scribe the interaction between the host immune system considering the role of 
humoral mediated immune responses, that is, the antibodies, and the cancer- 
ous cells. 

Let L, P, A and T be the number of large B-cells, plasma cells, antibodies and 
the cancerous cells respectively at any time t\. It has been observed that the 
number of antibody forming cells per spleen increased from the order of 10 2 
to 10 6 over 96 hours when a single injection of sheep red blood cells is given 
to mice [23]. However, when the absolute number of responding cells reaches 
the order of 10 7 , density dependent effects may prevent further growth of the 
lymphocyte population [24,25]. Therefore, we assume that the large B-cells 
grow logistically. The large B-cells can either proliferate with constant intrinsic 
growth rate a± and increase the B-cell population or they can undergo further 
differentiation into plasma cells at a constant rate b\. We assume a\ > b x to 
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insure that there can be a net growth in the large B-cells population. Thus, 
the governing equation for the large B-cells and the plasma cells are given by 



^ = a 1 «L(l--^)-6 1 (l-«)L, (1) 
dP 

— = 6 1 (l-«)L-|i 1 P, (2) 
ati 

where, u (0 < u < 1) is the fraction of the large B-cells which remains as 
the proliferating large B-cells and (1 - u) is the fraction that differentiate into 
plasma cells, K x is the carrying capacity of the large B-cell population and \X\ 
is the natural death rate of the plasma cells. 

It is clinically known that both the large B-cells as well as the plasma cells 
secrete antibodies, however, plasma cells secrete them at a much faster rate 
than the large B-cells. Therefore, the governing equation for the the antibody 
is given by 



= n L + r 2 P - fi 2 A - fcAT, (n < r 2 ), (3) 

(lt\ 

where r\ and r 2 are the rate at which the large B-cells and the plasma cells 
secrete antibodies respectively, \i 2 is the natural death rate of antibodies and 
f3i the rate at which the cancerous cells kill the antibodies. 

We assume that the cancerous cells grow logistically in the absence of the 
antibodies and the antibodies kill the cancerous cells directly [1]. Thus for the 
cancerous cells, the governing equation is 



f r rT(l-^)- k AT, (4) 

where r is the intrinsic growth rate, K 2 is the carrying capacity of the cancerous 
cells and fl 2 is the rate at which the antibodies kill the cancerous cells by direct 
interaction. The initial conditions for systems (1 — 4) are £(0) = Lq, P(0) = 
Po, A(0) = A and T(0) = T respectively. 



3 Parameter Estimation 



Appropriate parameter values determine the analysis and behavior of a math- 
ematical model to describe a given system. Therefore, to complete the de- 
velopment of our mathematical model, we estimate the values of the system 
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parameters in the following manner. The mean generation time for large B 
lymphocyte is approximated to be 6 to 48 hours [26]. Hence the growth rate 
ai of the large B lymphocyte is estimated to be between 0.02 and 0.2 hr' 1 [5]. 
Since the immune response to a T-independent antigen generally lacks any de- 
tectable immunological memory [27,28], the conversion of large lymphocytes 
back into small lymphocytes and the possible subsequent recycling of memory 
cells back into large lymphocytes have not been explicitly considered. The life 
time of plasma cells ranges from few days to few weeks. Thus the natural 
death rate /i 2 of the plasma cells is approximated to vary between 0.002 to 
0.02 hr- 1 [5]. 

The plasma and large B cells secrete antibodies at different rates. It is known 
that the plasma cells may be even 100 fold more active in immunoglobulin 
synthesis than the large B cells [29]. Thus, the rate at which the plasma cells 
secrete antibodies can be estimated to be 2 to 100 times that of the large B 
cells. The secretion rate of antibodies by a single cell also varies considerably. 
Different authors estimated the absolute rate of antibody secretion by a single 
cell. For example, Nossal and Makela [30] have been found in vitro that the rate 
at which a single cell secretes antibodies ranges from 100 to 1500 antibodies 
cell -1 sec.' 1 . Conrad and Ingraham [31] found in vivo values that range from 
8000 to 20000 antibodies cell' 1 seer 1 . 

The dynamics of cancerous cell (tumors) growth has been studied by numer- 
ous authors (V.A. Kuznetsov et. al [6], L.G. Pillis and A. Randunskaya [32], 
D. Krischner and J. C. Panetta [9], S. Banerjee and R.R. Sarkar [17]. They 
described the dynamics of tumor cell growth using the logistic growth function 
as dT/dt = rT(l — T/K2) in the absence of the immune responses, where r is 
the intrinsic growth rate and K 2 is the carrying capacity. L. G. Pillis et. al. 
[33] estimated the parameter values for r and K 2 by using the least-squares 
distance method and using optimization software with the data in Diefenbach 
et. al. [34]. In our case, we use r = 0.431 days -1 and K 2 = 9.8 x 10 8 cells 
from [33] for the analysis of our model. The interaction terms between anti- 
bodies and cancerous cells, namely, /3i and 2 and the natural death rate of 
antibodies, that is, ji 2 are estimated using synthetic data. 

To generate synthetic data, we consider the research article by L.G. de. Pillis 
et al. [33]. We first generate a figure (see fig. 2a) taking human data, patient 
9, as given in table 2 of [33]. This figure (see fig.2a) reflects the behavior of the 
model by De Pillis et. al. [33], who have used parameters taken from exper- 
imental results of patient from Rosenberg's study on metastatic melanoma. 
In the generated figure, they have investigated a 10 6 cells tumor, a tumor 
level that in silico innate immune system cannot control on its own. The fig- 
ure shows the effect of immunotherapy alone against the tumor, namely, a 
TIL (tumor infiltrating lymphocyte) injection, followed by short doses of IL- 
2. Using DataThief (www.datathief.org), the data for cancerous cells decay 
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are extracted (say, for 20 time points) and some noise is added to the data, 
which we name as observed values, T b s (say). To start the estimation process, 
we initially choose (within meaningful biological range) the values of the pa- 
rameters to be estimated, namely, /3 2 and /i 2 arbitrarily . Next, we solve 
the equations (1 — 4) numerically with these initial values of the parameters 
and obtain the solution of the model at those time points, where the observed 
values have been obtained, which we name as calculated values, T ca i (say). We 
now use least square method to minimize the sum of the residuals, namely, 
Yh1i(T^ — T^ s ) 2 to obtain the estimated values of the system parameters, 
namely, /3 2 and fi 2 . In practise, a MATLAB code has been developed to 
carry out the above process and the estimated values of the parameters are 
obtained as ft = 2.5448 x 10" 6 , ft = 2.4935 x 1(T 7 and fi 2 = 0.1277. Figure 
2b shows the best fit estimate for the model parameters. The corresponding 
estimated values of the parameters in non-dimensional form are respectively 
an = 5786.32, a 2 = 566.968 and r] 2 = 0.2963. 

We repeat these processes with five sets of data with different intensity of 
noises and succeeded in getting a range of values in which the parameters 
(3i, (3 2 and \i 2 lies. All the values of the system parameters are given in tabular 
form in Table 1. 



4 Analysis of the System 

To reduce the number of the system parameters and for numerical simulations, 
we non-dimensionalize the system using the following scaling, 



L P A T 
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r)2 = — , fci = —r^, k 2 = —r^i «i = , ol 2 

r rK 2 rK 2 r 
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Then the system can be described as 



dx 

~dl 

dy 

dt 

dz 

~dl 

dw 

~dl 



w{l — w) — a 2 zw, 



aux{l — x) — 6(1 — u)x, 



6(1 - u)x - r)iy, 



k\x + k 2 y — r] 2 z — a\zw, 



(8) 



(6) 



(7) 



(5) 
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with initial conditions x(0) = xo, y(0) = yo, z(0) = zq and w(0) = u>o- It can 
easily be shown that the right hand sides of system (5 — 8) are continuous and 
satisfy Lipschitz condition in Rl. 



4-1 Positivity and Boundedness of Solution 



Throughout this paper we denote qo := au — b{l — u) and assume that it is 
positive. 

Theorem 1 Let the initial conditions of system (5 — 8) be positive. Then the 
solutions (x (t), y (t),z (t), w (t)) of system are non-negative for all t > 0. 

Theorem 2 The solutions of system (5 — 8) with non negative initial condi- 
tions are bounded. 

Proof: From equations (5) and (6), we have 



^(x + y) = aux(l - x) - rny 

(ri 1 + au) 2 rj! + au 2 

= rj^x + y)- au{x 

Aau 2au 

(771 + au) 2 

Hence, using standard differential inequalities, we have x + y < • 

Similarly, from equations (7) and (8), since x(t) + y(t) is bounded (say by ko), 
we have 



— {z + w) = k\x + k 2 y — i] 2 z + w(l — w) — {a\ + a 2 )zw 

< k 2 {x + y) -r] 2 {z + w) - (w — ) + 

<k 2 k + (1 + ^ 2)2 - V2 (z + w). 



Hence (z + w) < + ■ Therefore, the solutions of system (5 — 8) are 

bounded. 
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4-2 Equilibria and Local Stability Analysis 



The equilibrium points of system (5—8) are E = (0, 0, 0, 0),Ei = (0, 0, 0, 1), E 2 = 

(x, y, z, 0) aud E* = (x, y, z*,w*), where x = y = z = 

and w* = 1 — a 2 z* such that z* is the positive roots of the quadratic equation 



aia 2 z 2 - (cci + r] 2 )z + (k x x + k 2 y) = 0. 



(9) 



Note: If q = au — 6(1 — u) < then the cancerous cells free E 2 and the 
positive interior E* equilibrium points of the system will not exist. But, from 
the linear stability analysis one can show that the equilibrium point E is 
unstable and the equilibrium point E\ is stable. 

We now study the linear stability of the cancerous cells free and the positive 
interior equilibrium points system (5 — 8). The jacobian matrix Je 2 at the 
boundary equilibrium point E 2 is given by 



J, 



i?2 



-(au-b(l-u)) 

6(1 -u) — r/i 

h k 2 -r] 2 -a x z 

1 - a 2 z J 

and the corresponding characteristic equation is 



(X + au- 6(1 - u))(X + 77i)(A + t? 2 )(A - (1 - a 2 z)) = 0. 



(10) 



Therefore, the equilibrium point E 2 (the cancerous cells free equilibrium point) 
is locally asymptotically stable (LAS) if 1 — a 2 z < 0. 

The Jacobian matrix J E * at the interior equilibrium point E* is given by 

\ 



J TP* 



(-(au-b(l -«)) 

6(1 - u) —77! 
ki k 2 —r] 2 — aiw* —a±z* 

y -a 2 w* -w* 



Then, the corresponding characteristic equation is 
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(X + au- 6(1 - u))(X + ^)[A 2 + (r] 3 + w*(l + e*i))A + (11) 

w*{j]2 + ol\ — 2a±a 2 z*)] = 

Theorem 3 Let z\ and z*_ be roots equation (6). The interior equilibrium 
point E*_ = (x,y, z*_,w*) , (the high number of cancerous cells equilibrium 
point) is LAS when ever it exists and = (x,y, z\,w*) is unstable. 

Proof: From the characteristic equation (11) the equilibrium point E* is LAS 
provided that the roots of the quadratic equation 



A 2 + (772 + w*(l + «i))A + w*(r] 2 + a x - 2a 1 a 2 z*) = 0, (12) 
are negative or with negative real parts. But from (9) we have, 



+ a-i — 2a 1 a 2 z* = =F\/ (^2 + «i) 2 - 4aia 2 (kix + k 2 y) 



Hence, the roots of equation (11) are negative or with negative real parts 
if i] 2 + «i — 2aia:2<2* > 0. Therefore, E*_ = (x,y, z*_,w*) is LAS and E* + = 
(x, y, z* +) w*) is unstable. 



4-3 Global Stability Analysis 



We now study the global stability of E 2 , the cancerous cells free equilibrium 
point and E*_, the high number of cancerous cells equilibrium point. 

Theorem 4 The cancerous cells free equilibrium point E 2 is globally asymp- 
totically stable ifa 2 > 4a %^lty) ' 

Proof: One can easily show that 



m («i + m)' 



k\x + k 2 y Aa\(k\x + k 2 y) 

Then, a 2 > 4^"^+^) implies that 1 — a 2 z < 0. Hence the cancerous cells 
free equilibrium point E 2 is LAS. From equation (5) we have 



x{t) = — - — ^° _ , where c=l — 



au(l — ce~ qot ) ' auxo 
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Hence, Hindoo x(t) = au ^ = x. From the definition of limit superior, for 
any e > 0, we have 



dy 

— = 6(1 - u)x - 7712/ < 6(1 - u)(x - e) - r]iy. 



Hence, 



6(1 — u)x 

hmsup?/(t) < = y, (since e is arbitrarily small). 

t->oo 1]i 

Similarly using the definition of limit inferior one can show that 



6(1 — u)x 
hmmt y[t) > = y. 

t^oo ww - ni 



Therefore, lim^ 00 y(t) = 6 ^ = y. On the other hand, from equation (7) 
and using the definition of limit superior, for e > 0, we have 



dz , , 

— = kix + k 2 y - rj 2 z - a x zw, 
at 



Therefore, lim sup z{t) < 



< kix + k 2 y - rj 2 z, 

, k\X + k 

>l2{ 

V2 

kix + k 2 y 



k x x + k 2 y 

<V2{ z). 

V2 



t^oo T) 2 



Thus, for every e > 0, there exist a time 7\ such that z(t) > (z — e) for t > 7\. 
Using this inequality in the equation (8) we have, 



— = w(l — a 2 z — w) < w(l — a 2 (z — e) — w) 

(JjL 

Since e is arbitrarily small and w(t) is non negative for all t > 0, we have, 
limsup^^ w(t) = 0. Therefore, for a 2 > ^"(l^+^y) ' ^ e cancerous ce U f ree 
equilibrium point E 2 = (x, y, z, 0) is globally asymptotically stable. 

Theorem 5 The positive interior equilibrium point E*_ is globally asymptot- 
ically stable if (1 — a 2 z) > 0. 
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Proof: For < a 2 < 4, the equilibrium point E*_ is LAS. From theorem (4), 
we have lim^oo x(t) = x, lim t ^ 00 ?/(t) = y and ^ < w(l — a 2 (z — e) — w). 
Then, limsup^^ w(t) < (1 — 012(2 — e)), if (1 — a 2 (z — e)) > 0. But since e > 
is arbitrarily small, we have lim sup t _ > . 00 w(t) < (1 — a 2 z) =: Wq, provided 
that (1 — a 2 z) > 0. From the definition of limit superior, for every e > 0, 
w(t) < (wq + e) for all t > 0. Substituting this in equation (7) we get 



dz 

— >k 1 x + k 2 y - (r) 2 + a-i(wl + e))z 
at 

= {m + «i(w + e))( ■ - z). 

Then, using differential inequality we have 

liminfzft) > l — ^- — — — -. 

Since, e > is arbitrarily small, we have 



v ■ t n\ ^ klX + ^ 
hmmt z(t) > =: z n . 

t^oo n 2 + aiWQ 

From the definition of limit inferior, for every e > 0, there exist a time T 2 such 
that z < Zq + e for all t >T 2 . Substituting this inequality in equation of (8) 
once again we get the inequality ^ > w(l — a 2 (zQ + e) — w), which intern 
implies that 



liminf w(t) > 1 — a 2 z^ =: w * 



t— >oo 



as e > is arbitrarily small. From the definition of limit inferior, for every 
e > 0, w > w* — e for all t > 0. Therefore, 



dz hx + k 2 y 

— <{V2 + a l {w -e)){—- — - -z) 

zt T) 2 + ai(w* — e) 



and hence 



kix + k 2 y 

hmsup z(t) < =: z , 

t^oo i] 2 + axw* 

since again as e > is so small. From the definition of limit superior again, for 
every e > 0, we have z > z* — e for all t > 0. Substituting this inequality in 
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equation (8), we get ^ < w(l — a 2 (z* — e) — w) and hence limsup^^ w(t) < 
1 — a 2 z* =: w*, as e is arbitrarily small. Now 



lim sup z(t) — lim inf z(t) = z* 



kix + k 2 y fax + k 2 y 

T] 2 + OtiW* T] 2 + atWQ 

a x {fax + k 2 y) 

(772 + aiw*)(772 + «i w o) 

a 2 ai(fax + fc 2 y) 
(7/2 + ai^oX^ + aiw*) 
a 2 a\(fax + k 2 y) 2 



V2(r]2 + aiWQ) 2 (r] 2 + aiw 



^0 - w ) 
2 o ~~ z o) 

-w*) < 0, 



which is a contradiction since Wq = 1 — > 0. Therefore, 
limsup^^ = liminf^oo z(t) = z*. Similarly, 

lim sup^oo w(t) —lim inf w(t) = Wq—w* = a 2 (zQ — z*) < 0, which is again 
a contradiction. Therefore, 



lim (x(t),y(t),z(t),w(t)) = {x, y, z*_ , w% 



and hence the positive interior equilibrium point E*_ = (x, y, z*_,w*) is globally 
asymptotically stable. 



4-4 Bifurcation Analysis 



In this section, we explore the critical parameter values where the qualita- 
tive behavior of the system changes. By using the Dulac-Bendixson theorem, 
one can show that system (5 — 8) has no closed orbit for positive solutions. 
Equations (5) and (6) can be evaluated analytically as, 

x(t) = — - — ^° _ , where c = 1 

au(l — ce qot ) auxo 

y (t) = y Q e- vlt + 6(1 - u) f e ni(s - t) x{s)ds, 

Jo 

which are not closed orbits. To show that the other solutions are also not a 
closed orbit, we consider the function m(w, z) = — . Then 

' \ ' / wz 
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d dz d dw 

9 1 9 1 

= tt( — ( k i x + k 2V - r) 2 z - ocizw)) + — ( — (w(l -w)- a 2 zw)) 
oz wz ow wz 

1 k 1 x + k 2 y 

z wz 



Since all the parameter values are positive, L < over the domain of interest 
and hence the system satisfies the Dulac-Bendixson theorem. Therefore, there 
are no limit cycles or homoclinic connections for the system. Similarly, from 
the values of the eigenvalues of the corresponding Jacobian matrix, there are 
no Hopf bifurcations, which may give rise to the occurrence of limit cycles. 

The system has different steady states depending on the values of the system 
parameters. The equilibrium points E = (0, 0, 0, 0) and E 1 = (0, 0, 0, 1) exist 
for all parameter values. The cancerous cells free steady state E 2 = (x, y, z, 0) 
exists if au — 6(1 — u) > 0. The other two positive steady states E* exist under 
the following conditions: 

(i) Let ai ~ V2 > 0. The equilibrium point E* + exist if \ < a 2 < ^"fc^+Cg) > 
and the equilibrium point E*_ exist if a 2 < 4 J"(l^l\ 2 ^ ■ 

(ii) Let «i — r] 2 < 0. The equilibrium point E*_ exists if a 2 < \. 
All these situations are sketched in fig. 3. 

The stability of the cancerous cells free equilibrium point E 2 changes as the 
value of the parameter a 2 passes through the critical value a 2 i = 4 = r^i+k^j ■ 
Hence a 2 , which represents the effectiveness of the antibodies to destroy the 
cancerous cells, is the bifurcation parameter for the system. The Jacobian 
matrix J at E 2 and its transpose have a simple eigenvalue A = with cor- 
responding eigenvector v T = (0, 0, —^f, 1) and w T = (0, 0, 0, 1) respectively 
at the bifurcation parameter value a 2 — 4. If f Q2 denote the vector of partial 
derivatives of the components of the right hand side of system (5 — 8) with 
respect to the scalar a 2 and Df(E 2 , a 2 )v is the directional derivative of f in 
the direction of v at the equilibrium point E 2 , then, 



w T i a2 {E 2 ) = 0, w T [Df a2 (E 2 , a 2 i)v] = = ™ and 

||v|| c^lMI 

Trn2p/P w m 2(0:1-772) 



w 1 [D%E ,a 2 i)(y,v) 



V2\\V\ 



which are non-zero for ai — r\ 2 7^ 0. Therefore, by Sotomayor theorem 
the system experiences transcritical bifurcation at the cancerous cells free 
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equilibrium point E 2 = (x,y,z,0), as the parameter a 2 passes through the 
bifurcation parameter value a 2 = a 2 i = 



Again, as the value of a 2 passes through the critical value a 2 = a 22 = 

4ai (fcit+fc 2 y) ' ^ e num t )er °f positive interior equilibrium points of the system 
changes. As the value of a 2 passes through the critical value a 22 , the pos- 
itive interior equilibrium points and E*_ collide and disappear. At the 
critical value of a 2 = a 22 the system has new interior equilibrium point E# = 
(x, y, 2a^al ' "lai 2 ) ' P rov ided a\—r}2 > 0. Then the Jacobian matrix J at E* and 
its transpose has a simple eigenvalue A = with corresponding eigenvectors 

V T = (0, 0, 1) and W T = (_ ^(&(l-u)fe 2 +mfci)(ai-^) _ a 2 fc 2 (a 1 - ??2 ) M _ 

ai ^ 2 ), 1) respectively, corresponding to the eigenvalue A = at the critical 
value a 2 = a 22 . Then, 



/ f a2 (^,a 22 ) = ^ 

w T [ J D 2 (f(^,a 22 )(v,v))] !?/: -' 



{Oil + V2)\\V\\ Z 

which are non-zero for a.\ 7^ r^ 2 . Therefore, by Sotomayor theorem, there is a 
smooth curve of the equilibrium point of the system inR xl passing through 
the point (E*, a 22 ) and tangent to the hyper plane 1R^_ x a 22 . Hence, the system 
has a saddle-node bifurcation as the parameter a 2 passes through the critical 
bifurcation parameter value a 2 = a 22 . 



5 Numerical Results and Biological Implications 



In this section we present the numerical results of system (5 — 8) for the 
parameter values given in Table 1. The steady sate E = (0, 0, 0, 0) and 
Ei = (0, 0, 0, 1) exist for the parameter values given in Table 1. The cancerous 
cells free equilibrium point E 2 = (x, y, z, 0) exists for the system parameters 
provided that au > 6(1 — u), which holds true for the parameter values given 
in Table 1. If a 21 = 0.227538 < a 2 < a 22 = 1111.032892, the positive interior 
equilibrium points E+, characterized by relatively low number of cancerous 
cells and E*_, characterized by relatively high number of cancerous cells exist 
keeping the other parameter values fixed. For < a 2 < a 2 \ = 0.227538, only 
the high number of cancerous cells equilibrium point E*_ exists. As the value 
of a 2 increases and passes through the critical value a 22 = 1111.032892, the 
equilibrium points E* + and E*_ collide and disappear. 

By linearizing the system about the steady states, the stability of the equi- 
librium points are determined, which is important from physiological point 
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of view. From the linear stability analysis we have seen that, the equilibrium 
points Eq = (0, 0, 0, 0) and E\ = (0, 0, 0, 1) are unstable for au — 6(1 — u) > 
and the cancerous cells free equilibrium point E 2 is stable if 1 — a 2 z < 0, that 
is, a 2 > \ = 0.227538 and unstable otherwise. The stability of E 2 implies that 
the antibodies can clear the cancerous cells. As a 2 passes through the criti- 
cal value a 2 \ = 0.227538, the stability of E 2 changes from unstable to stable 
and hence the system has a transcritical bifurcation at this point. Using the 
software Matcont, the point at which the system experiences transcritical 
bifurcation is BP = ( 0.550000 0.495000 4.394869 0.000000 0.227538 ), where 
x = 0.55, y = 0.495, z = 4.394869, w = 0.00 and a 2 = a 21 = 0.227538. As 
the value of a 2 passes through the critical value a 22 = 1111.032892 the two 
positive interior equilibrium points collide and disappear, where the system 
experiences a saddle node bifurcation at this point. Using Matcont, the point 
at which the system experiences the saddle node bifurcation is given by LP = 
( 0.550000 0.495000 0.000450 0.499974 1111.032892 ), where x=0.55, y=0.495, 
z=0.000450, w=0.499974 and a 2 = a 22 = 1111.032892. This is shown in fig.4a 
and 4b, where BP is the branching point, LP is the limit point and H is the 
neutral saddle, which has no meaning in the biological context. Also, it has 
been shown that E 2 attains global stability provided that a 2 > a 22 . 

The high number of cancerous cells equilibrium point E*_ is stable whenever it 
exists and is unstable always. In the region < a 2 < a 2 i, the high number 
of cancerous cells equilibrium point is the only stable (globally stable) equilib- 
rium point of the system and hence the cancerous cells succeed to survive. This 
is depicted in fig. 5, which implies that the cancerous cells will escape immune 
surveillance unless each and every cancerous cells are killed, that is, the system 
will ineluctably return to the high-cancer state if the treatment is stopped. 
Thus, in order to realistically effect a cure when the cancer free equilibrium 
point is unstable, a therapy or treatment must ensure that not only the can- 
cer burden must be reduced but the therapy itself is capable of changing the 
parameter values of the system. In this context, monoclonal antibody therapy 
of cancer [36] is suggested as treatment, which may be capable of changing 
the system parameters. Monoclonal antibodies (MAbs) can be used to target a 
number of cancer associated targets, including tumor associated blood vessels, 
vascular growth factors, diffuse malignant cells like leukemia, cancerous cells 
with a solid tumor and tumor associated stroma like fibroblasts. 

If a 2 \ < a 2 < a 22 , then both the cancerous cells free equilibrium E 2 and the 
high number of cancerous cells equilibrium E*_ are stable, that is, this region 
is the region of bistability. In this region, the system is sensitive to the initial 
conditions and some parameter values. Depending on the initial conditions 
and some of the values of the system parameters given in Table 1, a given 
initial value of the cancerous cell will either grow to the stable high number of 
cancerous cells equilibrium point E*_ (the cancerous cells succeed to survive) 
or decays to the cancerous cells free equilibrium point E 2 (the immune system 
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succeeds in eradicating the cancerous cells). Fig. 6 captures the dynamics of the 
sensitivity of the system with respect to the initial conditions. The sensitivity 
of the system on the parameter a 2 is shown in fig. 7. 

As the value of a 2 increases towards the critical value a 2 2, the quantity of the 
high number of cancerous cells equilibrium point decreases to about half the 
carrying capacity K 2 of the cancerous cells. When the value of a 2 increases 
and exceeds the critical value a 22) any amount of cancerous cells decays to 
zero after some time, that is, the immune system succeed in clearing the 
cancerous cells. The time at which a given amount of cancerous cells drops to 
zero depends on the initial conditions and some values of system parameters. 
When a 2 > a 22 , there is only cancer free equilibrium and it is stable. Fig. 8 
shows, once the critical point a 22 is crossed, the cancerous cells always decays 
to zero, no matter what the initial values are. 

As stated earlier, the effect of monoclonal therapy, which is capable of changing 
the system parameter, is evidently visible in this region, that is, when a 2 > a 22 . 
In this therapy, a monoclonal antibody can be directed to attach to certain 
cancerous cells. Cetuximab (Erbitux), a monoclonal antibody approved to 
treat colon cancer and head and neck cancers, attaches to receptors on cancer 
cells that accept a certain growth signal (epidermal growth factor). Cancer 
cells and some healthy cells rely on this signal to tell them to divide and 
multiply. Blocking this signal from reaching its target on the cancerous cells 
may slow or stop the cancer from growing. There are a number of monoclonal 
antibody drugs that are available to treat various types of cancer by (a) making 
the cancerous cells more visible to the immune system (b) delivering radiation 
to cancerous cells (c) slipping powerful drugs into cancerous cells. Recent 
clinical investigations show that they are capable of killing cancerous cells 
directly [1]. However, more clinical trials are necessary before it comes out in 
the form of medicine to mimic that effect. 



6 Conclusion 

In this paper, we have formulated a system of non-linear ordinary differential 
equations that describes the stimulatory effect of cancerous cells on immune 
cells in conjunction with antibodies. The model we have proposed is simple 
and more of general type. The major difference of our model from that of the 
existing literature in this direction is that the immune system we considered 
here is antibody mediated T-cell independent immune responses. In this dy- 
namics, a significant role is played by a 2 , the effectiveness of the antibodies 
to kill the cancerous cells directly. 

From our study we observe that for certain values of a 2 , one can control the 
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unlimited growth of the cancerous cells. From the analysis of our model, we 
determine some criterion for the existence of the equilibrium points of the 
system and the stability conditions. For a specific set of parameter values, 
we obtain five equilibrium points, namelythe zero equilibrium point E , the 
boundary equilibrium point Ei, the cancer-free equilibrium point E 2 and two 
positive interior equilibrium points E* . Initially, the cancerous cells free equi- 
librium point was unstable and the high number of cancerous cell equilibrium 
point was stable (a 2 < «2i)- This means any treatment to be effective, it 
must have the ability to change the system parameter and force this desirable 
equilibrium point to become stable. 

On the other hand, the stability of the high number of cancerous cells equi- 
librium point implies that reducing the cancer burden through any effective 
treatment is not sufficient enough to kill all the cancerous cells. Once the treat- 
ment stops, the system, even with an untraceable sign of cancer, will return 
to high number of cancerous cells state. However, alteration of system param- 
eters through monoclonal antibody therapy of cancer, may have the ability to 
change the stability nature of the cancer free equilibrium point and allowing 
a new treatment protocol to eradicate cancerous cells. In fact, it has been 
tested clinically that monoclonal antibodies directed to CD20 and HLA-DR 
can elicit homotypic adhesion followed by lysosome-mediated cell death in hu- 
man lymphoma and leukemia cells [36]. In the stability analysis of our model, 
we obtain a region of bi-stability and observe that the growth of the cancerous 
cells can be controlled and reduced at early stage of its growth by enhancing 
the patient's own defense mechanism, that is, the antibodies to fight against 
the cancerous cells. 

We perform bifurcation analysis of the system for the parameter a 2 , which 
represents the effectiveness of antibodies to eliminate cancerous cells. From 
the analysis, we observed that for certain parameter values of the system, the 
long term behavior of the system can be sensitive to the initial number of the 
cancerous cells. In the region where there is bi-stability (the stability region of 
both the cancer free and high number of cancerous cells equilibrium points), 
for the cancerous cell population that are very close to the boundary separat- 
ing the region of attraction of these equilibrium points, a slight change in the 
initial conditions will change the behavior of the system. One important impli- 
cation of this result for the treatment of cancer is that if both the cancer-free 
and high number of cancerous cells equilibria are stable, then for cancerous 
cell population that are close to the boundary separating the basins of attrac- 
tion of these equilibria, very small changes in the initial number of cancerous 
cells can have drastic consequences on the outcome of the cure. Therefore, by 
determining some of the values of the system parameters, we can achieve the 
stable cancerous cells free equilibrium point to cure the disease. At the end 
we would like to mention that the model can be modified further by adding 
delay terms, diffusive terms and stochastic fluctuations to the system to bring 



19 



out more rich dynamics of the system, which we propose as our future work. 
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Table 1 

Parameter values used for numerical illustration 



Parameters 


Values 


Scaled 


References 


cl\ (the growth rate of large B-cells) 


(0.02 - 0.2) hrr 1 


0.0464 -0.464 


f5l 


b\ (the conversion rate of large b-cells 
into plasma cells) 


0.01 hrr 1 


0.0232 


[5] 


/Xi(the natural death rate of plasma 
cells) 


(0.002 - 0.02) hrr 1 


0.00464-0.0464 


[51 


7^i(the carrying capacity of large B- 
cells) 


10 6 ce//s 




[51 


u (the fraction of daughter cells which 
remain as large B-cells) 


0.1 




[51 


ri(the rate at which large B-cells se- 
crete antibodies) 


100 Ab cell" 1 seer 1 


0.236754 


[29,30,31] 


r2(the rate at which plasma cells se- 
crete antibodies) 


1000 Ab cell" 1 seer 1 


2.36754 


[29,30,31] 


/^(the natural death rate of antibodies) 


0.1277 - 0.6465 seer 1 


0.2963 - 1.5 


estimated 


/?i(the death rate of antibodies due to 
interaction with cancerous cells) 


6.0436 x 10~ 9 - 2.5448 x 
10~ 6 cell^hrr 1 


13.7418 -5786.32 


estimated 


r (the intrinsic growth rate of cancerous 
cells) 


0.431 day- 1 




[33] 


K 2 (th.e carrying capacity of cancerous 
cells) 


9.8 x 10 8 cells 




[33] 


/^(the death rate of cancerous cells due 
to interaction with antibodies) 


9.0135 x 10~ 10 - 2.4935 x 
10~ 7 Ab' 1 hrr 1 


2.05 - 566.968 


estimated 
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Fig. 1. TTfce schematic diagram illustrates the interaction between large B-cells, 
plasma cells, antibodies and the cancerous cells 
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Fig. 2. The figure shows how syntectic data generated and used to estimate the 
parameters {5\ , /?2 and Hi . Fig. 2a is generated by taking the model and human data, 
patient 9 as given in table 2 of [33]. Fig. 2b shows the best fit estimate for the model 
parameters (3\ , /?2 and [12 ■ 
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Fig. 3. The figure shows the existence and stability regions of the equilibrium points 
E 2 , E\ andE*_, where k = kix+k 2 y. Ri (a%—r]2 < andO < a 2 < a 2 ± = 0.227538,) 
is the region where E 2 and E*_ exist such that E 2 is unstable and E*_ is sta- 
ble. R 2 (a± — rj2 > and < a 2 < a 2 \ = 0.227538), is the region where 
both E 2 and E*_ exist and stable (bistability region). R3 («i — rj 2 > and 
a 2 \ = 0.227538 < a 2 < a 22 = 1111.032892,) is the region where the equilibrium 
points E 2 ,E*_ and E\ exist and both E 2 and E_ are stable ( bistability region) and 
E+ is unstable. In the region R4 (a 2 > a 22 ) the equilibrium point E 2 exists and it 
is stable. 
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Fig. 4. The figure shows the points of bifurcation for the system parameter a,2, which 
has been obtained using the software Matcont. Here, H represents neutral saddle, BP 
is the branching point and LP is the limiting point. The branching point occurs at 
02 = 0.227538, which is identified as the transcritical bifurcation point and the limit 
point occurs at ai = 1111.032892 which is identified as the saddle node bifurcation 
of the system. Other system parameter values are fixed, which is given in Table 1. 
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Fig. 5. The figure shows that in the region < Q2 < «2i = 0.227538, the unique 
positive interior equilibrium point E*L characterized by high number of cancerous 
cells (nearly equal to the carrying capacity K%) is stable. Setting ai = 0.04 with any 
initial conditions (a) IC = (0.5,049,4.39,0.01) and (b)IC = (0.5,0.49,0.004,0.5), 
the cancerous cells grow to the high number of cancerous cells equilibrium point E*_ ■ 
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Fig. 6. The figure shows how the system is sensitive to the initial conditions in the 
region of bistability, that is, a<i\ < «2 < 022 • Setting 02 = 800, (a) initial conditions 
IC = (0.5,0.49,0.0004,0.3) shows that the cancerous cells grow to the high number 
of cancerous cells equilibrium point E*_. (b) However, a fraction of cell difference 
in the number of cancerous cells IC = (0.5,0.49,0.0004,0.2) shows that it decays to 
zero, to the cancerous free equilibrium point E%. 
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Fig. 7. T/ie figure illustrates that with the same initial condition IC = 
(0.5,0.49,0.0004,0.2) and for different values of 02, the system has different be- 
havior, (a) For ct<i = 726, the cancerous cells grow to the high number of cancerous 
cells equilibrium point and (b) for a<i = 727, it decay to the cancerous cells free 
equilibrium point. 
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Fig. 8. The figure shows when the value of «2 = 1200 greater than the criti- 
cal value «22 = 1111-03 for different values of the initial conditions (a) IC = 
(0.5, 0.49,0. 0004, 0.4)and (b) IC = (0.5,0.49,0.0004,0.85), the cancerous cells always 
decay to zero. 
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